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Abstract 

Network analysis has revealed several whole-network properties hypothesized to be general 
characteristics of ecosystems. These include pathway proliferation, and network non-locality, 
homogenization, amplification, mutualism, and synergism. Collectively these properties 
characterize the impact of indirect interactions among ecosystem elements. While ecosystem 
networks generally trace a thermodynamically conserved unit through the system, there appear 
to be several model classes. For example, trophic (TRO) networks are built around a food web, 
usually follow energy or carbon, and are the most abundant models in the literature. 
Biogeochemical cycling (BGC) networks trace nutrients like nitrogen or phosphorus and tend to 
have more aggregated nodes, less dissipation, and more recycling than TRO. We tested (1) the 
hypothesized generality of the properties in BGC networks and (2) that they tend to be more 
strongly expressed in BGC networks than in the TRO networks due to increased recycling. We 
compared the properties in 22 biogeochemical and 57 trophic ecosystem networks from the 
literature using enaR. We also evaluated the robustness of these results with an uncertainty 
analysis. The results generally support the hypotheses. First, five of the properties occurred in 
varying degrees in all 22 BGC models, while network mutualism occurred in 86% of the models. 
Further, these results were generally robust to a ±50% uncertainty in the model parameters. 
Second, the average network statistics for the six properties were statistically significantly 
greater in the BGC models than the TRO models. These results (1) confirm the general 
presence of these properties in ecosystem networks, (2) highlight the significance of different 
model types in determining property intensities, (3) reinforce the importance of recycling, and 
(4) provide a set of indicator benchmarks for future systems comparisons. Further, this work 
highlights how indirect effects distributed by network connectivity can transform 
whole-ecosystem functioning, and adds to the growing domain of network ecology. 

Keywords: network environ analysis; ecological network analysis; indirect effects; 
input-output analysis; systems ecology; network science 
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1 Introduction 

Complex systems are comprised of reticulated exchanges of energy, matter, and information 
that let members influence each other directly and indirectly. These connections can define 
member characteristics mmmm and create whole-system properties such as 
stability dEU and tolerance to failure [Hid- The degree of connectivity is important, but 
patterns of connectivity are critical for these effects UlMSIlj- Network models are ideal for 
mapping these exchanges, and network analyses characterize the subsequent relationships and 
quantify the connectivity patterns rnrnm- 

Ecosystems can be modeled as a network in which nodes represent species, functional groups 
of species, or non-living components, and directed edges trace the direct transactions of energy 
or matter among the nodes (Fig. [I}. Thus, this network model is like a road map for the 
physical exchange of resources. Ecologists then analyze the network model to understand the 
consequences of whole-ecosystem organization and the often hidden relationships among species 
and the external environment For example, Grami et al. [J5] found that 

in a food web model of Lake Pavin the fungal parasites of phytoplankton (chytrids) generally 
reduced the carbon lost from the pelagic zone. This effect was due in part because the chytrids 
increased the average path length that the carbon molecules traveled in the system, which they 
further hypothesized created greater system stability. Small et al. [91] used ecosystem networks 
to compare nitrogen cycling in the Laurentian Great Lakes, and discovered that Lake Superior 
had a higher denitrification efficiency than Lake Huran or Lake Erie despite its having a low 
areal denitrification rate. This ecosystem network analysis (ENA) builds our understanding of 
the importance of connectivity patterns in ecosystems, contributes to the rise of network 
ecology Ena, and more broadly contributes to the emerging domain of network science and its 
approach to studying complex systems in general [7l [22if70l[72] . 

While ecosystem networks are formed by tracing a thermodynamically conservative unit 
through the system [22] , there appear to be several classes of ecosystem network models with 
distinctive features mmm- Trophic networks (TRO) are the most common class of ecosystem 
networks. They are built around a food web, and typically trace either energy or carbon. The 
South Carolina oyster reef model [9] and Chesapeake Bay ecosystem models [7] are typical 
examples of this class. A second class of ecosystem network model focuses on biogeochemical 
cycling (BGC), and typically traces a key nutrient like nitrogen or phosphorus. Christian et 
al. [24] suggested that by comparison to the trophic models, biogeochemical cycling models tend 
to have more aggregated nodes, less dissipation and more recycling. Models of nitrogen cycling 
in the Laurentian Great Lakes EH, Neuse River Estuary [4], and the Cape Fear River 
Estuary PE] are examples of this class. In this work, we investigate the generality of selected 
ecosystem network properties in the two most abundant classes (TRO and BGC). 

1.1 Six Hypothesized Characteristics of Ecosystem Networks 

Through the development and application of ENA, several hypotheses regarding the 
organization of whole ecosystems have emerged JISElEillEE! • Here, we focus on six of these 
hypotheses: pathway proliferation, network non-locality, network homogenization, network 
amplification, network mutualism, and network synergism. 

Pathway proliferation occurs when the number of pathways in an ecosystem tends to rapidly 
increase as path length increases Emms]. In this literature a pathway is defined as a 
sequence of nodes and edges (e.g., Detritus —► Detritivore —> Consumer in Fig. [TJi), and the 
length of the pathway is the number of edges m, where to = 2 in the example. This is 
sometimes called a walk in graph theory as both nodes and edges can be repeated [20]. These 
longer pathways create a reticulated web of possible routes of influence. Pathway proliferation 
is illustrated in Fig. [if’. 

Network non-locality is the tendency for indirect flows (flows over pathways where to > 1) to 
cumulatively exceed direct flows in ecosystems @3ii75][ES] ■ This occurs despite the diminishing 
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Figure 1: Hypothetical ecosystem network model, (a) The directed graph represents an 
ecosystem with the three nodes (Detritus, Detritivores, and Consumer) along with the internal 
flows fij and boundary inputs Zi and outputs yj of energy or matter. The matrices (b) summarize 
the steady-state dynamics of energy or matter fluxes [M L' 2 T' 1 ] in the system if we assume a 
10% trophic transfer efficiency. Pathway proliferation (c) is shown as the number of pathways 
(walks) from node 2 to 1 ( 012 ) and from node 3 to 2 ( 032 ) increase with path length, and (d) 
shows the expected rapid decline in flow intensity () between nodes as path length increases. 


flow intensity as path length increases (Fig. |TJl) due to the thermodynamically open nature of 
ecosystems. The implications of this property are that organisms not only have the ability to 
influence each other indirectly, but that these indirect relations may be both quantitatively and 
qualitatively dominant. Further research has shown that these indirect interactions have the 
power to transform the nature of the relationships among species. For example, Bondavalli and 
Ulanowicz m found that due to indirect interactions across the trophic network in the 
Everglades the American alligator was more beneficial to some of its prey like frogs than they 
were detrimental. This complicates environmental management and restoration jl01l74| . 

Network homogenization is the tendency of the indirect flows to more evenly distribute 
resources in the system. Patten et al. m observed that when they applied ENA flow analysis, 
the integral (boundary + direct + indirect) flow intensities were more evenly distributed than 
the direct flow intensities alone. They inferred that the indirect flow intensities must be 
spreading the resource exchange. 

Network amplification is defined as “the increase of energy-matter ultimately utilized at 
destination compartments compared to the quantities originally introduced at source 
compartments” [50] ■ Though it may at first seem counter intuitive, this property is made 
possible by energy and matter recycling J5D1EI] ■ It implies that some exchanges are quite 
important for the system function. 

Network mutualism and Network synergism both reflect the tendency for indirect 
interactions to make the relationships among the network nodes more positive than they appear 
from the direct interactions [35l[39l[82], as with the previously mentioned case of the frog and 
alligator. Network mutualism examines the number of positive and negative interactions, while 
network synergism compares the magnitudes of the positive and negative interactions. 
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1.2 Existing Evidence 

These hypothesized general ecosystem properties were first identified in relatively small, 
well-connected, trophic-based ecosystem models (e.g., oyster reef model 0), but subsequent 
research has shown that four of the six properties appear to be general in trophic models. For 
example, pathway proliferation is known to occur in all well connected models m- Higashi and 
Patten [48j[49] showed algebraically why network non-locality is likely to occur in ecosystem 
networks. Fath and Patten |40j introduced a statistic to quantify network homogenization and 
showed that the property occurred in two ecosystem models and a set of arbitrary networks. 
Subsequently, Fath [34] and Fath and Killian [38j found evidence to support the network 
non-locality and homogenization hypotheses in large-scale hypothetical model ecosystems built 
from an ecosystem assembly algorithm; however, they found that network amplification 
occurred only rarely. Salas and Borrett |86j found that the network non-locality hypothesis held 
in 74% of 50 empirically-based trophic ecosystem networks, and Borrett and Salas [18] showed 
that network homogenization occurred in 100% of the same 50 models. Fath and Patten |39] 
provided an algebraic proof for why network synergism should always occur in ecosystem 
models. Thus, the existing evidence supports the claim that pathway proliferation, network 
non-locality, network homogenization, and synergism are general properties of trophic-based 
ecosystem models, while network amplification is an occasional property or possibly a modeling 
artifact. Network mutualism has rarely been investigated. 

The majority of evidence for these ecosystem properties is based on trophic ecosystem 
models. We posit that each ecosystem model class may have different characteristic network 
properties or intensities of the properties. That said, the ENA framework predicts that all six 
described ecosystem properties are general across the model classes because these models still 
represent an ecosystem at their core. 

1.3 Objectives 

We investigated whether these six ENA properties generally occur in biogeochemical cycling 
ecosystem networks. We specifically tested two hypotheses. First, we hypothesized that the 
pathway proliferation, non-locality, homogenization, amplification, mutualism, and synergism 
properties tend to occur in the BGC models because the properties are purportedly general to 
ecosystems. Second, we hypothesized that these six properties would tend to have a greater 
intensity in BGC models than in TRO models because (1) BGC models tend to have higher 
rates of recycling, and (2) all of the properties tend to increase with recycling [SUlZZliSS] • If 
these hypotheses hold, then this work extends the generality of the hypothesized ecosystem 
properties across both the trophic and biogeochemical cycling model classes, and provides 
evidence that differences between model classes can be functionally significant. 

2 Materials and Methods 

2.1 Ecosystem Network Selection 

Our approach to testing the generality of the ecosystem properties was to apply ENA to a set of 
ecosystem network models and compare the results. As such, this study has characteristics of a 
systematic review and met a-analysis, and we can characterize our model selection for the study 
with an adaptation of the PRISMA (Preferred Reporting Items for Systematic reviews and 
Meta-Analyses) guidelines [57]. Figure [5] provides an overview of this decision process. 

To identify ecosystem networks for this study, we started with two existing collections. The 
first set of 48 models was compiled by Dr. Robert Ulanowicz 

http://www.cbl.umces.edU/~ulan/ntwk/network.h.tml over the course of his career 
developing and promoting ENA. The second is a collection of 100 ecosystem models extracted 
from the literature and distributed with the enaR R package for ENA M- These were collected 
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Figure 2: Information flow for ecosystem network model selection adapted from the PRISMA 
(Preferred Reporting Items for Systematic reviews and Meta-Analyses) guidelines. 


from colleagues and periodic literature searches (2000 - 2014). Together, this identified 124 
unique ecosystem networks. In an initial screening step we excluded two of these networks 
because the original model sources were unavailable for review. 

We then assessed the eligibility of the remaining 122 ecosystem networks for our study. To 
be included, the models must have met five inclusion criteria. First, the models needed to fit 
the characteristics of a trophic or biogeochemical cycling network. We excluded two models 
based on this criteria - one was modeled a network of water flow in a watershed, and another 
was for the Polish economy. Second, we included only the highest resolution (least aggregated) 
networks for a particular system. This criteria is necessary because some systems have been 
modeled multiple times with differing network order (number of nodes used to partition the 
system). Third, to be included in our study, the model had to be empirically-based in the sense 
that the investigators were attempting to model a specific ecosystem and a portion of the model 
was parameterized with empirical measurements from that system. Fourth, for analytical 
reasons (see next section) the models included had to be at steady-state (inputs equal outputs). 
Thus, we excluded models that were not initially at steady-state that could not be 
automatically balanced using enaR. Fifth, there were several cases where there were multiple 
models of the same system for different times or conditions. If an averaged model was also 
initially identified for the site, we included only the averaged model to limit the bias of any 
particular system. In some cases a similar conceptual model was parameterized for different 
elements (e.g., C, N, P, Mg). We included this variety because we expected this to provide 
alternative views of the same system. 

Through this process we selected 79 ecosystem network models from 41 records. This 
includes 22 biogeochemical cycling ecosystem networks (Table SI) and 57 trophic networks 
(Table S2). Eight distinct systems are modeled in the set of biogeochemical network class and 
there are 40 separate systems in the trophic network class. This set of trophic models overlaps 
(49 of 50) with those analyzed for network non-locality and network homogenization (18 ) 186] . 
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2.2 Ecological Network Analysis 

While the specifics of the ENA algorithms have been presented in the literature j25 (141ll51H89 _ 
and made available in software packages mmm, the techniques are not well known. Further, 
in some cases we modified the traditional ENA network statistics to facilitate their 
consideration across networks with different numbers of nodes. Thus, this section describes the 
analyses used in detail. 

2.2.1 Ecosystem Network Models 

Like network models of any complex system, ecosystem networks represent the system with a 
set of nodes and edges (Fig. [1]). Nodes represent biological species, groups of species, or 
non-living components. For example, a living node could represent detritivores, which is a 
functional grouping of species based on their eating habits, and non-living nodes might be 
particulate organic matter (POM) or inorganic nitrogen in the form of ammonium (NH 4 + ). 
Edges are directed and weighted. They indicate the flow of energy or matter from one node j to 
another i over some time period, F = [fij\ i,j = 1, 2,..., n. These networks allow only one edge 
from any node to any other node, and self-loops may occur (e.g., cannibalism in a food web). 
These models typically focus on tracing a single thermodynamically conservative tracer. Thus, 
an edge may capture the effect of many physical, chemical, and biological processes including 
feeding, excretion, death, and chemical conversions. 

Ecosystems are thermodynamically open systems. This implies that the systems have inputs 
and outputs from their environment, which are modeled as a second class of directed edges. 
Input edges, z = [zi] do not have a specified starting node and output edges y = [yj\ have no 
terminal node. Alternatively we could assign a node to represent the environment, but this is 
uncommon in the ENA literature. 

Given an ecosystem network, we can apply ENA. There are several distinct analyses within 
ENA pmifRU : we applied structural, flow, and utility analyses. 

2.2.2 Structural Analysis 

Structural analysis in ENA focuses on properties of the network topology. This topology is 
specified by the existence of connections among the nodes; node and edge weights are ignored. 
The analysis also typically focuses on the within system connections, so the boundary inputs 
and outputs are also ignored. Thus, the analysis primarily operates on the binary adjacency 
matrix, A = [a^], induced from F. If there is a flow from j to i, aij = 1; otherwise, a. t j = 0. 

A number of network statistics are used to characterize the digraph structure. The most 
common are the total number of nodes n (graph order), the number of edges or links 
L = and connectance C = L/n 2 (network density). These statistics are a necessary 

starting point, but they do not describe the pattern of connections among the nodes. 

Structural analysis is also concerned with features of the network connectivity, which are 
characterized with graph theory and matrix algebra [HE], a strongly connected network is one 
in which every node is reachable by a pathway of some length from every other node. Some 
networks can be decomposed into one or more strongly connected components (SCC). These 
SCCs are the largest subset of nodes that are still reachable from each other [9]. Most food 
webs and ecosystem networks are composed of one or sometimes two SCCs with more than one 
node, which we call large SCCs (T2MT9]. Often these networks also have a set of additional 
nodes that are weakly connected to the others; they are only reachable if we ignore the edge 
directions. For example, a plant node in a trophic network tracing carbon typically has an edge 
showing consumption by one or more herbivores; however, there is typically no input of carbon 
to the plants from within the system since they take it up from the atmosphere. We can 
quantify these aspects of network connectedness by counting the number of large (n > 1) SCCs 
in the network and the percent of network nodes that participate in the SCCs. 
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Pathway proliferation appears when we count the number of pathways or walks between any 
two nodes as we increment the pathway length m (Fig.QJi) [76lf80| . Raising the adjacency 
matrix to the power m returns a matrix A m whose elements a^ 1 ' 1 provide the walk count from j 
to i. In strongly connected networks, a( m ' ) —> oo as m —> oo. 

The rate of pathway proliferation varies among ecological networks mm- This variation is 
due to the size of the network, its direct connectance, and the topological arrangement of the 
direct edges Ej. The asymptotic rate of pathway proliferation in an SCC is equivalent to the 
dominant eigenvalue of the adjacency matrix associated with the SCC, Ai(A) jT2]. If a network 
contains only one SCC, then the pathway proliferation rate for the network and the SCC are 
identical. Ai(A) may also serve as a measure of the topological cyclicity in the network [321136] . 
As Ai(A) scales between 0 and n, we used a re-scaled pathway proliferation rate 
PPR = (Ai(A))/n to compare the pathway proliferation rates among networks of different 
sizes. Notice that 0 < PPR < 1. 

2.2.3 Flow Analysis 

Flow analysis characterizes the geometry of the ecological networks, taking into consideration 
the flow weights [4611551 . Further, it requires the boundary flux information. Given its 
origination in economic Input-Output Analysis, flow analysis can have either an input or output 
orientation; we focus on the output orientation for this paper. The input analysis typically 
generates qualitatively similar results m 

The ENA flow analyses we applied have four key steps. The first step is to calculate the 
node throughflows, which are a sum of the total amount of energy or matter flowing into or out 
of each node. This can be calculated as 


T m put = ^ + Zi ^ j = i j 2,..., n ), and (1) 

3 

n 

^output = £ fij + y. (f , j = 1, 2,..., n). (2) 


The flow analysis applied here assumes the network model is at steady-state, which implies that 

^-linput _ rji Output _ rj~i 

i j 3 

The second step in flow analysis is to calculate the direct flow intensities among nodes 
(G = \gij])- For the output orientation, these values are defined as 


9ij — 


Tj 


(3) 


Elements of G indicate the flow intensity from j to i over pathways of length rn = 1. 

The third step in our analysis is to find the integral flow intensity matrix N = [n^]. This 
matrix satisfies the following equation. 


T = Nz 


(4) 


The riij elements are effectively the sum of the flow intensities from j to i across pathways of all 
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Table 1: Formulation of six network statistics to characterize hypothesized general properties of 
ecosystems. 


ENA Property 

Symbol 

Statistic 

Pathway Proliferation 

PPR 

Ai(A)/n 

Network Non-locality 

I/D 

£((N-I-G)z)/£Gz 

Network Homogenization 

HMG 

CV(G)/CV{N) 

Network Amplification 

AMP 

> 1) / {n{n — 1)) 

Network Mutualism 

MUT 

(#«ij > 0) / (#Vij < 0) 

Network Synergism 

SYN 

(£ £ v ij e ® >0 ) / ( £ £ Vij e R<o) 


lengths (m = 0,1,2, 3,, oo). Thus, we can re-represent N as the sum of the following series: 

OO 

N = G m (5) 

m—0 

= + G 2 + ... + G m + ... . (6) 

Boundary Direct Indirect 

The elements of G m indicate the flow intensity from j to i over all pathways of length in. 

Notice that when equation [6] is substituted back into equation U the elements of G°z map the 
boundary inputs into T, G 1 ^ are the direct flows, and G m z where m > 1 are the indirect flows. 
Given that ecosystems are thermodynamically open and the specification of our model, this 
infinite series is convergent. Therefore, we can find N using the identity (I — G) , where I is 
the identity matrix. 

From this primary flow analysis, we derive several whole-network statistics relevant to this 
paper. The first two are based on the throughflow vector. Total system throughflow (TST) is 
the sum of the node throughflows (TST = ^G T)) and is an indicator of the size or activity of 
the system [U3[55]- Average path length ( APL = TST/Y^Zi) indicates the average amount of 
activity generated by a unit of input into the system |43j . This is similar in concept to the 
multiplier effect in economics |88| and has been termed network aggradation in ecology [511] • 

The third indicator we used is the Finn Cycling Index (FCI), which indicates the fraction of 
throughflow derived from recycling as opposed to acyclic flows pnHFPT| . 

The remaining flow-based network statistics quantify three of the remaining ENA properties 
of interest (Table [0. Network non-locality is characterized by the ratio of indirect-to-direct 
flows in the network (//£>); when I/D > 1, indirect effects are said to dominate direct effects 
and network non-locality is present. Network homogenization is a ratio of the coefficient of 
variation in the direct flow intensity matrix to the coefficient of variation in the integral flow 
matrix ( HMG ). When HMG > 1 the flows are more evenly distributed (lower variance) in the 
integral than in the flow matrix, which implies that the indirect flows are more evenly 
distributing the system activity across the network and network homogenization is present. 
Network amplification occurs when an off-diagonal element in the integral flow matrix exceeds 
1, which indicates that the flow is amplified from the initial boundary input. The strength of 
amplification ( AMP ) is quantified by counting the number of non-diagonal elements are greater 
than one. To allow comparison across network models of different sizes, we divide this count by 
the total number of possible off diagonal elements. 
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2.2.4 Utility Analysis 

The network mutualism and network synergism properties appear from the application of ENA 
utility analysis. This technique determines the net relationships among nodes in the network 
that are derived from the direct transactions ISSMHiiM]- Here, we use the flow based utility 
analysis [JUST]- The first step in this analysis is to calculate the direct utility matrix as 

D = [dij] = fi \ fj \ i,j = 1,2, ...,n (7) 

-L i 

such that —1 < dij < 1. Further, d,, :] shows the direct utility intensity from i to j over pathways 
of length m = 1. As in flow analysis, the integral utility intensity matrix is then computed as 

OO 

U = ^2 D m (8) 

m—0 

= + JD^+ p 2 + ... + D m + ... . (9) 

Boundary Direct Indirect 

The elements of U indicate the integral utility intensity from i to j. When the dominant 
eigenvalue of D is less than one, this power series converges such that U = (I — D) _1 . When 
the convergence criteria fails, the identity relationship should not be mathematically 
appropriate. Despite this mathematical problem, use of the identity has been common in the 
literature and it appears to return ecologically meaningful results. Why this appears to work 
despite the mathematical issue is an open problem in ENA. As such, in this paper we complete 
the utility analysis even when the convergence criteria fails, but identify the results as 
mathematically suspect. 

The next step in utility analysis is to re-scale and dimensionalize both the direct and 
integral utility matrix by multiplying them by the node throughflows, 

A = [8ij] = TD, and (10) 

T = K-] = TU, (11) 

where T is a n x n matrix with the node throughflows on the principle diagonal and zeros 
elsewhere. Thus, Sij is the dimensionalized [M L' 2 T -1 ] direct utility i receives from j, and Vij is 
the dimensionalized integral utility. 

The tests for network mutualism and network synergism involve the ratios of positive and 
negative utilities in U, which have been termed the benefit-cost ratios (Tabled]). In both cases, 
the property is said to occur if the ratios are greater than 1, indicating that the positive utilities 
exceed the negative utilities. 

2.3 Data Analysis 

Our primary data analysis for this paper had three main components. First we used the enaR 
package to apply ENA to the selected BGC models to test the hypothesized generality of the 
properties in the biogeochemical cycling ecosystem networks. We evaluated the presence of each 
ecosystem property as expressed in the network statistics described (Table |T]). Second, we used 
a two-sample Wilcoxan signed rank test to compare the intensity of each ecosystem property 
between the BGC and TRO sets of models (significance criterion a = 0.05). Third, to test the 
effect of model size and recycling on the ENA properties, we applied the non-parametric 
Spearman Rank correlation. 
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2.4 Uncertainty Analysis 

To evaluate the effect of flux uncertainty on the ENA results, we applied an inverse linear 
modeling technique EDIISS1- We used the limSolve R package [H] to generate an ensemble of 
10,000 plausible models for each network by simultaneously but independently perturbing the 
model fluxes by a random amount up to ±50% drawn from a uniform distribution. These 
perturbed models were constrained to meet the steady-state assumption of ENA. We then 
applied ENA to the ensemble of models and calculated the 95% confidence interval for each 
network statistic. 


3 Results 

The results of this study generally support the hypotheses. We report the results in four stages. 
We first present the results of the ENA network statistics in the biogeochemical cycling models, 
and then we show the values for the trophic based models. The third subsection compares the 
distribution of these properties between the two model classes. In part four, we confirm the 
importance of cycling intensity for these six ecosystem properties. 

3.1 ENA properties in Biogeochemical Cycling Ecosystem Networks 

Pathway proliferation, non-locality, homogenization, amplification, and synergism occurred to 
varying degrees in all 22 BGC models analyzed (Fig. [3]). Network mutualisms occurred 86% of 
the models. Thus, the results support the first hypothesis that these are general properties of 
ecosystems. 

The relative rate of pathway proliferation ranged from 0.11 to 0.47 (Fig. Eh) - This 
proliferation rate reflects the results for the whole network because the models each contain a 
single large (n > 1) strongly connected component. Furthermore, all of the nodes participate in 
the model’s large SCC except in the Hubbard Brook Ca model in which only 3/4 do. As this is 
purely a topological property, the flux uncertainty analysis did not affect these results. 

In all of the BGC, models I/D > 1 (logio(I/D) > 0). This indicates the dominance of 
indirect flows and the presence of network non-locality. The median realized I/D was 12.2 in 
the BGC models, but it ranged between 2.07 in the polyhaline model of N cycling in the Cape 
Fear River Estuary to 201.1 in the Mn cycling model in the Puerto Rican rain forest model. 

The 95% Cl intervals from the uncertainty analysis (error bars) indicate a range of uncertainties 
for the specific values for each network, but none of these intervals cross the interpretation 
threshold of I/D = 1. 

Likewise, the homogenization statistic exceeded one in all BGC models ( HMG > 1). This 
indicates that the integral flow intensity matrices were less variable than the direct flow 
intensity matrices. When the pathways of increasing length were considered, the particular 
nutrient being traced through the network was more evenly distributed between the nodes. 
HMG ranged from 1.84 to 13. Again, the 95% Cl do not cross the interpretation threshold of 
HMG = 1. 

Network amplification occurred (AMP > 0) in all 23 BGC models; the median value faction 
of possible off-diagonal positions with elements greater than 1 was 0.75. The lowest level of 
amplification was 2% in the polyhaline Cape Fear Estuary model of the N cycle and 
amplification peaked at 100% in nine of the networks. The 95% CIs reinforced this result. 

The network mutualism statistic ranged from 0.77 to 4 in the BGC models, where 86% of 
them had a value larger than one. Only the Chesapeake Bay N model and the Sylt-Rpmp Bight 
N and P models failed to exhibit the network mutualism property. Notice, however, that the 
95% uncertainty interval for the Chesapeake Bay N model includes one. As expected all of the 
models exhibited network synergism. The synergism statistic ranged from 0.55 to 1.78 on a log 
scale. However, the Utility matrix failed the convergence criterion (Ai(D) < 1) in all but three 
models, making these results mathematically suspect (Fig. [3]) . 


10/05] 
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Figure 3: Intensity of the (a) relative pathway proliferation rate ( PPR) (b) network non-locality 
(I/D), (c) network homogenization (HAIG), (d) network amplification (AMP), (e) network 
mutualism (MUT), and (f) network synergism (SYN) in 22 biogeochemical cycling ecosystem 
networks. Notice that non-locality and synergism were log 10 transformed to facilitate comparisons 
(recall that log 10 (l) = 0). Models are ordered from smallest (left) to largest (right) n, and the 
error bars represent the 95% confidence interval based on the uncertainty analysis. The open 
bars in network mutualism and network synergism indicate that the utility analysis failed to meet 
the convergence criteria. 
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3.2 ENA Properties in Trophic Ecosystem Networks 

There is also general support for five of the six ENA hypotheses in the 57 TRO models (Fig. 0}, 
though these results are tempered by the uncertainty analysis. Note that the uncertainty 
analysis technique failed to build an ensemble of perturbed models for 14 of the TRO models 
due to “incompatible constraints” that made the solutions impossible to identify |2S], including 
the eight largest models. 

Pathway proliferation occurred in all but two models: the Silver Springs and the English 
Channel. The maximum relative pathway proliferation rate was 0.4. Non-locality occurred in 
75.4% of the initial models, including all 4 of the models not previously considered. However, 
with the 50% flux uncertainty the 95% Cl of the non-locality statistics for 14 of the models 
crossed the existence threshold, which leaves us uncertain of their status. Network 
homogenization occurred in all 57 TRO models with values ranging from 1.04 to 1.96. Network 
amplification was found in 75% of the trophic models, with values ranging from 0.01% to 17%. 
Network synergism occurred in all 57 models with logio transformed parameter values ranging 
from 0.38 to 1.02. The 95% Cl for network homogenization, amplification, and synergism all 
exceeded the exceeded the property existence thresholds. 

The evidence for the generality of network mutualism was not strong; the network statistic 
ranged from 0.60 to 3 in the TRO models. The network mutualism parameter only exceeded 
one, indicating the ENA property presence, in 28 of the 57 models (49%). Further, the 
uncertainty analysis showed that the evidence was inconclusive for 7 of these 28 models. Also 
notice that 47.4% of the TRO models failed to meet the convergence criterion for the utility 
power series, which makes these results mathematically suspect. 

The results for non-locality and homogenization are not surprising given the known results 
on 49 of these models [181186]. However, these two properties were found in all seven of the 
unique models in this set. 

3.3 BGC versus TRO network models 

When we compared the distributions of the ENA statistics between the two classes of ecosystem 
models (BGC vs TRO), we found a significant difference for each metric (Fig. [5])• As expected 
from our class definitions, these distributions indicate that the BGC models in our data set 
tended to have fewer nodes, connectance, and a higher degree of recycling. They also tended to 
have a higher average path length. Further, as we hypothesized the relative pathway 
proliferation rate, network non-locality, network homogenization, network amplification, 
network mutualism, and network synergism tended to be greater in the BGC models. A 
non-parametric Wilcoxan rank sum test showed a statistically significant difference between 
BGC and TRO models for the ten metrics compared (Table [2]). The p-values are all below the 
critical value of 0.05, which indicates that there is sufficient evidence to reject the statistical 
null hypotheses of no difference. 

3.4 Importance of Cycling 

We hypothesized that the explanation for why the ENA properties would be more intense in the 
BGC models was due to their tendency to have greater cycling. Figure [6] confirms for this data 
set that all six properties tend to increase with the Finn Cycling Index, though the shape and 
intensity of the relationship varies. The correlations between FCI and the ENA property 
indicators was statistically significant (Table [3|. 
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Figure 5: Comparison of nine network properties in a set of biogeochemical cycling (BGC, N = 22) 
and trophic (TRO, N = 57) ecosystem networks. These properties include (a) number of nodes, 
(b) pathway proliferation rate, (c) network non-locality, (d) Finn Cycling Index, (e) network 
homogenization, (f) network amplification, (g) average path length, (h) network mutualism, and 
(i) network synergism. 


4 Discussion 

4.1 Ecological Significance 

Ecosystem network analysis and its associated environ theory have been developing since the 
mid 1970s [TTli4Tl|46l[83] . Today, these techniques and theory appear to be in a maturation 
phase. This work makes two primary contributions to this maturing systems ecology 
understanding of ecosystems. 

First, the results presented here extend the evidence in support of five of the hypothesized 
general ecosystem properties: pathway proliferation, network non-locality, network 
homogenization, network amplification, and network synergism. Theses properties tend to 
occur in empirically-based network models of two classes: trophic and biogeochemical cycling. 
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Table 2: Summary of Wilcox Rank Sum statistical comparison of the distribution of ten network 
statistics between BGC and TRO models. 


Network Statistic 

W 

p-value 

Significance 

n 

240 

1.28e-05 

* 

c 

912 

3.92e-03 

* 

FCI 

1,218 

6.05e-10 

* 

APL 

1,118 

3.22e-07 

* 

PPR 

985 

2.35e-04 

* 

I/D 

1,243 

1.06e-10 

* 

HMG 

1,260 

3.09e-ll 

* 

AMP 

1,234 

1.72e-10 

* 

MUT 

1,041 

5.86e-06 

* 

SYN 

1,163 

4.72e-06 

* 


Table 3: Results of a Spearman rank correlation test between the Finn Cycling Indices and six 
Ecosystem Network Analysis statistics for 79 ecosystem models. 


Network Statistic 

rho 

t 

p-value 

Significance 

PPR 

0.51 

37588.97 

2.69e-06 

* 

I/D 

0.94 

5211.53 

8.40e-37 

* 

HMG 

0.74 

21543.63 

8.77e-15 

* 

AMP 

0.93 

5811.33 

4.85e-35 

* 

MUT 

0.62 

31595.33 

1.58e-09 

* 

SYN 

0.71 

23824.64 

2.36e-13 

* 


Network mutualism was more variable, especially in the BGC models, and strength of network 
amplification did decline with network size as it did in simulated ecosystem networks [34] . This 
evidence adds to the work from studies of empirically-based ecosystem networks [33|7j3[|86] as 
well as generative ecosystem models (SUES]. Further, these results were generally (but not 
always) robust to the uncertainty introduced into the flux estimates. 

Second, the results suggest that the lens used to view the ecosystem - trophic or 
biogeochemistry - can influence the intensity of the observed network properties. In general, 
models with more cycling tend to have a more intense expression of these properties. This 
concern with analytical point of view is similar to Reiners’ [86] recognition of multiple 
complimentary models for ecosystems. 

When considered in aggregate, these six hypothesized general properties of ecosystems reveal 
some of the ecological consequences of ecosystem connectivity through the lens of a single 
transactive currency. The intensity and pattern of direct connections among ecosystem 
members leads to indirect interactions that alter the ultimate nature of the interactions. 

4.2 Limitations and Opportunities 

A limitation of this study concerns the two sets of ecosystem network models. While the models 
themselves likely vary in quality, there are fewer BGC models in our sample than TRO models 
and the BGC models tended to have fewer nodes. This likely effects our inference about the 
differences between the two network classes; however, we expect the general trends reported in 
this paper to hold. This claim is supported by the general increase in the ecosystem network 
properties in the nitrogen and phosphorus models when compared to the related carbon models 
of both the Chesapeake Bay and Sylt-Rpmp Bight models (Figs EE]). 

To determine the robustness of our results to potential parameter uncertainty, we applied an 
uncertainty analysis based on inverse linear modeling [60][92l[99]. While this kind of global 
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Figure 6: Scatter plots showing the relationship between (a) relative pathway proliferation rate, 
(b) network non-locality (I/D), (c) network homogenization (HMG), and (d) network amplifica¬ 
tion (AMP) in 79 ecosystem network models. 


uncertainty analysis [87] is useful for estimating our confidence in the results, it has some 
limitations. First, the technique assumes the model parameters are independent when we might 
expect that some of the parameters may have natural inter-dependencies. Because the 
technique misses the potential covariances, the technique as applied is conservative, leading to 
potential type II statistical inference errors. Second, in the absence of better parameter 
uncertainties, we applied a uniform ±50% parameter perturbation. In reality, some of the 
parameters may be more and less certain [7]. The effect of this heterogeneity of natural 
variability and uncertainty is an issue that needs further exploration. 
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While this second issue could increase or decrease the estimated 95% confidence intervals, 
we expect that when combined with the first issue that the confidence intervals are over 
estimates of the true uncertainty. If true, this bolsters the evidence presented in this paper for 
the support of the six ecosystem properties. A third issue that will affect the potential broad 
applicability of the inverse linear modeling based uncertainty analysis is the occasional failure of 
the procedure to identify models that meet the mathematical constraints. This is a known 
challenge for the technique [99] . 

Despite the methodological challenges, uncertainty analyses are essential for making 
informed inferences from ENA ®mm- Thus, we expect that uncertainty analyses of some 
form will become standard practice for future ENA work. 

4.3 Network Science 

Network models and analyses are useful in a wide range of disciplines including medicine and 
genetics [3011105] . ecology [2l U5ll29] . sociology [100] . education [69] . and business [28ll53lll03j . 
From these domain specific theoretical and methodological developments, a network science is 
emerging j22[[ZQl[ 72.. Brandes et al. [22] define network science broadly as the “study of the 
collection, management, analysis, interpretation, and presentation of relational data” as used to 
study a wide range of complex systems. While some expect that this approach will identify 
properties, structuring forces, and organizational consequences of connectance common to 
complex systems in general [7], Brandes et al. [22] claim a more modest goal for network theory. 
They suggest that network theory will focus on how to abstract systems into useful network 
representations and suggest appropriate ways of querying these models. 

Ecosystems have all of the hallmarks of complex adaptive systems [65[[53]. This paper 
highlights a network model to investigate ecosystem energy and matter transactions and what 
appear to be common properties of this system class. We suspect that this has a broader 
relevance for network science. This abstraction and subsequent analysis was initially built upon 
economic Input-Output Analysis j46j|64], but has developed substantially to meet the demands 
of ecological and environmental questions. Thus, this ecosystem network model and analysis 
maybe relevant for other kinds of complex systems that can be modeled similarly. 

The network science type proposition that ecosystem network models can be divided into 
two (or more) distinct classes is pivotal for the work presented in this manuscript. The evidence 
in Fig. [5] suggests that the models generally follow the expected differences between the 
proposed classes. Further, the hypothesized systems properties do seem to be general in both 
groups. Other possible classes include hydrologic models m and urban metabolism or 
industrial material flow models mmmm- 

Unfortunately, the distinction between TRO and BGC models is not always crisp. For 
example, the core conceptual model for the Sylt-Rpmp Bight TRO ecosystem was only slightly 
modified to build the BGC models. For example, the node representing dissolved organic 
carbon (DOC) was re-purposed to represent dissolved nitrogen, and then the edge weights were 
estimated for the nitrogen fluxes rather than the carbon fluxes. Thus, the model topology is 
still largely that of the original trophic model, while the geometry responds to the different 
currency. In many cases this kind of transition makes sense because multiple nutrients 
(currencies) travel together in some ecosystem transactions. For example, when deer eat grass 
they obtain energy, carbon, and nutrients like nitrogen and phosphorus through the same 
process, though the stoichiometric ratios may vary in time and space [621193] . Thus, some of the 
TRO and BGC models are not strictly independent, and perhaps in the future they could be 
combined through the multilayer network framework recently introduced [60] . Despite this 
issue, we assert that the ecosystem model classification is a useful analytical technique that can 
capture important differences among the ecosystem model types. We look forward to the 
extension of these ideas to other classes of ecosystem models like the urban metabolism 
models [23][59][6^, or other kinds of input-output systems [6511101] . 
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Table SI: Twenty-two empirically-based ecosystem network models of biogeochemcial cycling. 


Model 

Units 

n 

C 

Input 

TST 

FCI 

Reference 

Hubbard Brook (Ca)(Waide) 

kg Ca Ha -1 yr — 1 

4 

0.25 

11 

168 

0.76 

El 

Hardwood Forest, NH (Ca) 

kg Ca Ha - ^ yr — ^ 

4 

0.31 

11 

200 

0.80 

0 

Duglas Fir Forest, WA (Ca) 

kg Ca Ha - ^ yr — ^ 

4 

0.31 

4 

54 

0.74 

0 

Duglas Fir Forest, WA (K) 

kg K Ha - ^ yr 1 

4 

0.31 

0 

45 

0.97 

0 

Puerto Rican Rain Forest (Ca) 

kg Ca Ha - ^ yr — ^ 

4 

0.31 

43 

274 

0.57 

0 

Puerto Rican Rain Forest (K) 

kg K Ha -1 yr -1 

4 

0.31 

20 

433 

0.86 

0 

Puerto Rican Rain Forest (Mg) 

kg Mg Ha -1 yr — 1 

4 

0.31 

10 

70 

0.58 

0 

Puerto Rican Rain Forest (Cu) 

kg Cu Ha" 1 yr” 1 

4 

0.31 

0 

2 

0.37 

0 

Puerto Rican Rain Forest (Fe) 

kg Fe Ha - ^ yr — ^ 

4 

0.31 

0 

7 

0.95 

0 

Puerto Rican Rain Forest (Mn) 

kg Mn Ha -1 yr -1 

4 

0.38 

0 

7 

0.98 

0 

Puerto Rican Rain Forest (Na) 

kg Na Ha - ^ yr — ^ 

4 

0.31 

64 

140 

0.24 

m 

Puerto Rican Rain Forest (Sr) 

kg Sr Ha” 1 yr” 1 

4 

0.31 

0 

1 

0.71 

0 

Tropical Rain Forest (N) 

g N m- 2 d- 1 

5 

0.24 

10 

71 

0.48 

0 

Neuse River Estuary (N, AVG) 

mmol N m — 2 season - ^ 

7 

0.45 

795 

41,517 

0.89 

0 

Cape Fear River Estuary (N, oligohaline) 

nmol N cm - ^ d — ^ 

8 

0.36 

3,802 

7,088 

0.20 

a 

Cape Fear River Estuary (N, polyhaline) 

nmol N cm - ^ d — ^ 

8 

0.36 

3,068 

5,322 

0.17 

m 

Lake Lanier (P) Averaged 

mg Pm -2 day -1 

11 

0.21 

95 

749 

0.40 

0 

Baltic Sea (N) 

mg N m — ^ day - ^ 

16 

0.15 

2,348 

44,510 

0.67 

0 

Chesapeake Bay (N) 

mg N m -2 yr -1 

36 

0.12 

73,430 

484,325 

0.33 

2 ' 

Chesapeake Bay (P) 

mg P m -2 yr -1 

36 

0.12 

9,402 

101,092 

0.51 

ITo] 

Sylt-Rpmp Bight (N) 

mg N m -2 yr -1 

59 

0.09 

99,613 

363,693 

0.23 

m 

Sylt-R 0 m 0 Bight (P) 

mg P m -2 yr -1 

59 

0.09 

2,510 

57,708 

0.66 

m 
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Table S2: 57 empirically-based, trophic ecosystem networks. 


Models 

Units 

n 

C 

Input 

TST 

FCI 

Reference 

Marine Coprophagy (oyster) 

kcal m~^ yr —1 

4 

0.25 

379 

549 

0.12 

po! 

Lake Findley 

gC m“ 2 yr- 1 

4 

0.38 

21 

50 

0.30 

1201 

Mirror Lake 

gC m“ 2 yr" 1 

5 

0.36 

72 

217 

0.32 

I20l 

Lake Wingra 

gC m 2 yr 1 

5 

0.40 

478 

1,517 

0.40 

1201 

Marion Lake 

gC m — 2 yr -1 

5 

0.36 

87 

242 

0.31 

1201 

Cone Springs 

kcal m -2 yr - 

5 

0.32 

11,819 

30,626 

0.09 

1241 

Silver Springs 

kcal m — 2 yr - 

5 

0.28 

21,296 

29,175 

0.00 

1161 

English Channel 

kcal m — 2 yr - 

6 

0.25 

1,096 

2,280 

0.00 

m 

Oyster Reef 

kcal m —2 yr — ^ 

6 

0.33 

41 

83 

0.11 

GD 

Baie de Somme 

mgC m —2 d — 1 

9 

0.30 

876 

2,034 

0.14 

I 2 TI 

Bothnian Bay 

gC m -2 yr -1 

12 

0.22 

44 

183 

0.23 

1221 

Bothnian Sea 

gC m- 2 yr- 1 

12 

0.24 

117 

562 

0.31 

1221 

Ythan Estuary 

gC m“ 2 yr" 1 

13 

0.23 

1,258 

4,181 

0.24 

Si 

Sundarban Mangrove (virgin) 

kcal m -2 yr - 

14 

0.22 

117,959 

441,214 

0.16 

P9l 

Sundarban Mangrove (reclaimed) 

kcal m — 2 yr — ^ 

14 

0.22 

38,484 

103,056 

0.05 

1191 

Baltic Sea 

mg C m —2 d — 1 

15 

0.17 

603 

1,973 

0.13 

m 

Ems Estuary 

mg C m- 2 d _ 1 

15 

0.19 

282 

1,067 

0.32 

m 

Southern Benguela Upwelling 

mg C m —2 d 1 

16 

0.23 

714 

2,545 

0.31 

- 

Peruvian Upwelling 

mg C m —2 d 1 

16 

0.22 

14,927 

33,491 

0.04 

m 

Crystal River (control) 

mg C m- 2 d _ 1 

21 

0.19 

7,357 

15,062 

0.07 

[25l 

Crystal River (thermal) 

mg C m —2 d — 1 

21 

0.14 

6,018 

12,032 

0.09 


Charca de Maspalomas Lagoon 

mg C m —2 d — 1 

21 

0.12 

1,486,230 

6,010,331 

0.18 

m 

Northern Benguela Upwelling 

mg C m~ 2 d _1 

24 

0.21 

2,281 

6,608 

0.05 

m 

Swartkops Estuary 

mg C m~ 2 d _ 1 

25 

0.17 

2,859 

8,949 

0.27 

rm 

Sunday Estuary 

mg C m~ 2 d _ 1 

25 

0.16 

4,441 

11,939 

0.22 

1231 

Kromme Estuary 

mg C m~ 2 d 1 

25 

0.16 

2,571 

11,087 

0.38 

1231 

Okefenokee Swamp 

g dw m~ 2 y 1 

26 

0.20 

2,533 

12,855 

0.48 

1301 

Neuse Estuary (early summer 1997) 

mg C m — 2 d — 1 

30 

0.09 

4,385 

13,827 

0.12 

m 

Neuse Estuary (late summer 1997) 

mg C m —2 d — 1 

30 

0.11 

4,639 

13,035 

0.13 

3: 

Neuse Estuary (early summer 1998) 

mg C m ~ 2 d — 1 

30 

0.09 

4,568 

14,025 

0.12 

m 

Neuse Estuary (late summer 1998) 

mg C m —2 d 1 

30 

0.10 

5,641 

15,031 

0.11 

m 

Gulf of Maine 

g ww m~~ 2 yr -1 

31 

0.35 

5,053 

18,381 

0.15 

P 21 

Georges Bank 

g ww m -2 yr -1 

31 

0.35 

4,380 

16,889 

0.18 

1121 

Middle Atlantic Bight 

g ww m -2 yr -1 

32 

0.37 

4,869 

17,916 

0.18 

1121 

Narragansett Bay 

mgC m -2 yr -1 

32 

0.15 

693,845 

3,917,246 

0.51 

1151 

Southern New England Bight 

g ww m~ 2 yr -1 

33 

0.35 

4,717 

17,597 

0.16 

1121 

Chesapeake Bay 

mg C m -2 yr -1 

36 

0.09 

888,791 

3,227,453 

0.19 

7 

Mondego Estuary (Zostera sp. Meadows) 

g AFDW m~~ 2 yr -1 

43 

0.19 

4,030 

6,822 

0.03 

17I18| 

St. Marks Seagrass, site 1 (Jan.) 

mg C m —2 d — 1 

51 

0.08 

514 

1,315 

0.13 

a 

St. Marks Seagrass, site 1 (Feb.) 

mg C m —2 d — 1 

51 

0.08 

601 

1,590 

0.11 

m 

St. Marks Seagrass, site 2 (Jan.) 

mg C m- 2 d _ 1 

51 

0.07 

602 

1,383 

0.09 

m 

St. Marks Seagrass, site 2 (Feb.) 

mg C m~ 2 d _ 1 

51 

0.08 

800 

1,921 

0.08 

s 

St. Marks Seagrass, site 3 (Jan.) 

mg C m -2 d _ 1 

51 

0.05 

7,809 

12,651 

0.01 

0 

St. Marks Seagrass, site 4 (Feb.) 

mg C m~ 2 d — 

51 

0.08 

1,432 

2,865 

0.04 

0 

Sylt-Rpmp Bight 

mg C m —2 d 1 

59 

0.08 

683,448 

1,781,028 

0.09 

m 

Graminoids (wet) 

g C m- 2 yr- 1 

66 

0.18 

6,272 

13,676 

0.02 

[551 

Graminoids (dry) 

g c m — 2 yr -1 

66 

0.18 

3,472 

7,519 

0.04 

1291 

Cypress (wet) 

g C m — 2 yr -1 

68 

0.12 

1,418 

2,571 

0.04 

1261 

Cypress (dry) 

g C m~ 2 yr -1 

68 

0.12 

1,035 

1,919 

0.04 

1261 

Lake Oneida (pre-ZM) 

g C m~ 2 yr -1 

74 

0.22 

1,034 

1,697 

0.00 

1131 

Lake Oneida (post-ZM) 

g C m — 2 yr -1 

76 

0.22 

810 

1,462 

0.00 

1131 

Bay of Quinte (pre-ZM) 

g C m“ 2 yr- 1 

74 

0.21 

989 

1,517 

0.00 

1141 

Bay of Quinte (post-ZM) 

g C m — 2 yr -1 

80 

0.21 

1,163 

2,107 

0.01 

1141 

Mangroves (wet) 

g C m — 2 yr -1 

94 

0.15 

1,531 

3,265 

0.10 

1281 

Mangroves (dry) 

g C m~ 2 yr -1 

94 

0.15 

1,531 

3,272 

0.10 

128') 

Florida Bay (wet) 

mg C m — 2 yr — ^ 

125 

0.12 

738 

2,720 

0.14 

1271 

Florida Bay (dry) 

mg C m — 2 yr — ^ 

125 

0.13 

547 

1,778 

0.08 

1271 
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